C Pacific sardine simulation code. Written by F. Hurtado-Ferro and A.E. Punt
C Please reference the authors if you use this code.
C
      IMPLICIT NONE
      INCLUDE "SARDINE.INC"
C
C     Local variables 
      INTEGER Year,Isim
C      REAL*8 TAC
C      
C     Read in the specifications
      CALL ReadInBiol()
C
C     Read in the parameters WRIof the control rule

      OPEN(UNIT=14,FILE="SUMMARY.OUT")

C      ISEED1 = -89910
C      ISEED2 = -76901
C      ISEED3 = -99701
C      ISEED4 = -65244	  
	  
      DO 10000 Isim = 1,Nsim
C
C      Initialize
       CALL Initialize(Isim)
C
       DO 11000 Year = 1,Nyear
C        TAC = 0
        CALL Proj1Yr(Isim,Year)
11000  CONTINUE        
C
10000 CONTINUE     
      CALL SUMMARY()
      CLOSE(14)
C
C      WRITE(*,*) Ninit
C
      STOP
      END
C
C ===============================================================================================
C
      SUBROUTINE ReadInBiol()
C
      IMPLICIT NONE
      INCLUDE "SARDINE.INC"
C
C     Local variables
      INTEGER Age
C
      OPEN(UNIT=13,FILE="Sardine.Dat")
      READ(13,*)
      READ(13,*) Nage
      WRITE(*,'(A5,1x,I3)') "Nage",Nage
      READ(13,*)
      READ(13,*) (Weight90s(Age),Age=0,Nage)
      READ(13,*)
      READ(13,*) (Weight80s(Age),Age=0,Nage)
      READ(13,*)
      READ(13,*) (Matu(Age),Age=0,Nage)
      READ(13,*)
      READ(13,*) (Fecu(Age),Age=0,Nage)
      READ(13,*)
      READ(13,*) (Ninit(Age),Age=0,Nage)
      READ(13,*)
      READ(13,*) (Selex(Age),Age=0,Nage)
      READ(13,*)
      READ(13,*) (SelxMC1(Age),Age=0,Nage)
      READ(13,*)
      READ(13,*) (SelxMC2(Age),Age=0,Nage)
      READ(13,*)
      READ(13,*) (SelxPNW(Age),Age=0,Nage)
      READ(13,*)
      READ(13,*) Natm
      READ(13,*)
      READ(13,*) AlogRPS, BlogRPS, PhilogRPS
      READ(13,*)
      READ(13,*) Alpha, Beta1, Beta2
      READ(13,*)
      READ(13,*) EpsInit,ProwR,SigmaR,Rerr,Rtype
      READ(13,*)
      READ(13,*) Vinit,Gtype,UseG,Verr
      READ(13,*)
      READ(13,*) Ybar,Period,ProwV,SigmaV,Phi,MeanSST, Amplitude
      READ(13,*)
      READ(13,*) SigmaB, Berr
      READ(13,*)
      READ(13,*) Chi
      READ(13,*)
      READ(13,*) Power
      READ(13,*)
      READ(13,*) SigmaI,Ierr
      CLOSE(13)
C
      OPEN(UNIT=13,FILE="Sardine.Spec")
      READ(13,*)
      READ(13,*) Nyear
      READ(13,*)
      READ(13,*) Nsim
      READ(13,*) 
      READ(13,*) Emsymin, Emsymax	  
      READ(13,*) 
      READ(13,*) FracMin, FracMax
      READ(13,*) 
      READ(13,*) Cutoff
      READ(13,*) 
      READ(13,*) Maxcatch
      READ(13,*) 
      READ(13,*) OFLopt
      READ(13,*) 
      READ(13,*) HCRopt
      READ(13,*) 
      READ(13,*) CutOpt
      READ(13,*) 
      READ(13,*) ISEED1, ISEED2, ISEED3, ISEED4
      READ(13,*) 
      READ(13,*) TVselx
      READ(13,*) 
      READ(13,*) TVweight
      READ(13,*) 
      READ(13,*) MultFleet
      READ(13,*) 
      READ(13,*) HyperSt, HPSa
      CLOSE(13)
C      
      RETURN
      END
C
C ===============================================================================================
C
      SUBROUTINE Initialize(Isim)
C
C This subroutine sets up variables which are common across a run
C
      IMPLICIT NONE
      INCLUDE "SARDINE.INC"
C
C     Global variables
      INTEGER Isim
C
C     Local variables
      INTEGER Iage,Year
      REAL*8 NN(0:MaxAge),S0,TwoPi,EpsV,EpsR,Term1,Term2
      REAL*8 SSBV,Bio1P,Idev
      REAL*8 XNORM
      REAL*8 ANcount1, ANcount2
      REAL*8 Qweight
      EXTERNAL XNORM
C
      TwoPi = 6.283185
C
C     Initial N-vector
      ANcount1 = 0      
      ANcount2 = 0
      AgeN(Isim,0) = 0
      DO 10000 Iage = 0,Nage
       N(0,Iage) = Ninit(Iage)
       ANcount1 = ANcount1 + N(0,Iage)*Iage
       ANcount2 = ANcount2 + N(0,Iage)
10000 CONTINUE
      AgeN(Isim,0) = ANcount1/ANcount2
      RecRep(Isim,0) = N(0,0)       
C
C     By default, set the weight vector equal to that of 1991-2000
      DO 14000 Iage = 0,Nage
       Weight(Iage) = Weight90s(Iage)
14000 CONTINUE
C
C     Compute the SSB and 1+ biomass
      SSBV = 0
      Bio1P = 0
      DO 11000 Iage = 1,Nage
C       SSBV = SSBV + Fecu(Iage)*Matu(Iage)*N(0,Iage)
       Bio1P = Bio1P + Weight(Iage)*N(0,Iage)
11000 CONTINUE
      DO 12000 Iage = 2,Nage
       SSBV = SSBV + Weight(Iage)*N(0,Iage)
12000 CONTINUE
      SSB(Isim,0) = SSBV
      B1p(Isim,0) = Bio1P
C
C     Generate the environmental variables
      Vdev(0) = Vinit
      V(0) = Vdev(0) + MeanSST
      Vrep(Isim,0) = V(0)
C	  
      DO 31000 Year = 1,Nyear
C
C      Create G
       Term1 = Sin(TwoPi*(Year-Ybar)/Period)
       IF (MOD(Year,NINT(Period)).LT.Period/2.0) THEN
        Term2 = 1
        ELSE
         Term2 = -1
       ENDIF
C       WRITE(*,*) Year,Term1,Term2
       IF (Gtype.EQ.1) THEN
        G(Year) = Term1
        ELSE
         G(Year) = Term2
       ENDIF
C
C      Update
       EpsV = XNORM(1,0.0d0,1.0d0,ISEED1)*SigmaV
       EpsV = EpsV * Verr
       Vdev(Year) = ProwV*Vdev(Year-1)+(1-ProwV)*G(Year)*Amplitude
       Vdev(Year) = Vdev(Year) + SQRT(1-ProwV**2)*EpsV
       V(Year) = Vdev(Year)*UseG + MeanSST
       Vrep(Isim,Year) = V(Year)
C       V(Year) = MeanSST
C       WRITE(*,*) Year,G(year),V(Year)
C       WRITE(*,*) Vdev(Year)
C
31000 CONTINUE      
C
C     Generate the rec_dev
      Eps(0) = EpsInit
      DO 32000 Year = 1,Nyear
       EpsR = XNORM(2,0.0d0,1.0d0,ISEED2)*SigmaR
       EpsR = EpsR * Rerr
       Eps(Year) = prowR*Eps(Year-1)+SQRT(1-ProwR**2.0)*EpsR
C       WRITE(*,600) Year,"E",EpsR,Eps(Year),V(year)
32000 CONTINUE
C      
C     Add measurement error to the environmental variable
      DO 33000 Year = 0,Nyear
       Idev = XNORM(1,0.0d0,1.0d0,ISEED3)*SigmaI*Ierr
       I(Year) = V(Year)+Idev
C       WRITE(*,'(A3,1x,F8.5,1x,A3,1x,F8.5,1x,A4,1x,F8.5)') "V",V(Year),
C     +       "Iv",I(Year),"Idev",Idev
33000 CONTINUE
C 
C     Add measurement error to the 1+ biomass (i.e. an "assessment")
      EpsB(0) = 1
      DO 34000 Year = 1,NYear   
       EpsB(Year) = XNORM(1,0.0d0,1.0d0,ISEED4)*SigmaB - 
     +                  (SigmaB**2)/2.0
C       WRITE(*,*) EpsB(Year),Berr 
       EpsB(Year) = EXP(EpsB(Year)*Berr)
34000 CONTINUE	   
C     
C     SENSITIVITIES
C     1. Time varying weight at age
      IF (TVweight.EQ.1) THEN
C      Compute Q
       Qweight = -33.245934 + 2.115014*V(0)
       Qweight = MIN(1.0,Qweight)
       Qweight = MAX(0.0,Qweight)
C      Compute the SSB and 1+ biomass
       DO 50000 Iage = 0,Nage
        Weight(Iage) = Qweight*Weight80s(Iage)+
     +                 (1-Qweight)*Weight90s(Iage)
50000  CONTINUE
       SSBV = 0
       Bio1P = 0
       DO 51000 Iage = 1,Nage
C        SSBV = SSBV + Fecu(Iage)*Matu(Iage)*N(0,Iage)
        Bio1P = Bio1P + Weight(Iage)*N(0,Iage)
51000  CONTINUE
       DO 52000 Iage = 2,Nage
        SSBV = SSBV + Weight(Iage)*N(0,Iage)
52000  CONTINUE
       SSB(Isim,0) = SSBV
       B1p(Isim,0) = Bio1P
      ENDIF
C
C     For debugging only
C      DO 50000 Year = 1,Nyear
C      WRITE(*,600) Year,G(Year),V(Year),I(Year),Eps(Year)
C50000 CONTINUE       
C
      RETURN
600   FORMAT(1x,I4,1x,A1,1x,7(F7.3,1x))      
      END
C
C ===============================================================================================
C
      SUBROUTINE Proj1Yr(Isim,Iyr)
C
C This subroutine updates the basic dynamics
C
      IMPLICIT NONE
      INCLUDE "SARDINE.INC"
C
C     Global variables
      INTEGER Iyr,Isim
      REAL*8 HG
      REAL*8 HarvGuide	  
      REAL*8 Eff, ObsEff
C
C     Local variables
      INTEGER Age,II
      REAL*8 FF,Z,Fmin,Fmax,Cnum,Cpred
      REAL*8 SSBV,Bio1P,RecR,SSBflat
      REAL*8 ANcount1, ANcount2
      REAL*8 ACcount1, ACcount2
      REAL*8 Lselx,Jselx,Qweight,MaxSelx
      REAL*8 qind
      REAL*8 RickerA, RickerB
      RickerA = -4.598
      RickerB = -3.457e-07
      IF(Rtype.EQ.4) THEN
       RickerA = -3.6978
       RickerB = -1.4680e-06	
      ENDIF	   
C      
C      WRITE(*,600) Iyr,(N(Iyr-1,Age),Age=0,Nage)
C
C     SENSITIVITIES
C     1. Time varying weight at age
      IF (TVweight.EQ.1) THEN
C      Compute Q
       Qweight = -33.245934 + 2.115014*V(Iyr)
       Qweight = MIN(1.0,Qweight)
       Qweight = MAX(0.0,Qweight)
C      Compute the weight-at-age vector
       DO 50000 Age = 0,Nage
        Weight(Age) = Qweight*Weight80s(Age)+
     +                (1-Qweight)*Weight90s(Age)
50000  CONTINUE
      ENDIF
C     For debugging
C      WRITE(*,'(I3,1x,F5.3,1x,20(F5.3,1x))') Iyr,Qweight,
C     +                                       (Weight(Age),Age=0,Nage)
C      WRITE(*,'(I3,1x,F6.3,1x,F5.3)') Iyr,V(Iyr),Qweight
C
C     SENSITIVITIES
C     2. Time varying selectivities
      IF (TVselx.EQ.1) THEN
C      Compute L
       Lselx = -31.407868 + 1.998082*V(Iyr)
       Lselx = MIN(1.0,Lselx)
       Lselx = MAX(0.0,Lselx)
C      Compute J
       Jselx = -23.790353 + 1.518223*V(Iyr)
       Jselx = MIN(1.0,Jselx)
       Jselx = MAX(0.0,Jselx)
C      Compute the MexCal selectivity for the year
       DO 51000 Age = 0,Nage
        SelxMCy(Age) = Lselx*SelxMC1(Age) + 
     +                 (1-Lselx)*SelxMC2(Age)
51000  CONTINUE
C      Normalize (set the max selx to 1)
       MaxSelx = MAXVAL(SelxMCy)
       DO 52000 Age = 0,Nage
        SelxMCy(Age) = SelxMCy(Age)/MaxSelx
52000  CONTINUE
C      Now, compute the total selectivity
       DO 53000 Age = 0,Nage
        Selex(Age) = Jselx*SelxMCy(Age) +
     +               (1-Jselx)*SelxPNW(Age)
53000  CONTINUE
C      And normalize again
       MaxSelx = MAXVAL(Selex)
       DO 54000 Age = 0,Nage
        Selex(Age) = Selex(Age)/MaxSelx
54000  CONTINUE
      ENDIF
C
C     Do the "Assessment"
      Bhat = B1p(Isim,Iyr-1)*EpsB(Iyr-1)
C	  
C     SENSITIVITIES 
C      Hyperstability
      IF (HyperSt.EQ.1.0) THEN
        qind = 1
        qind = HPSa*1/SQRT(B1p(Isim,Iyr-1))
        qind = MAX(qind,1.0)
        Bhat = qind*B1p(Isim,Iyr-1)*EpsB(Iyr-1)
C        WRITE(*,'(I3,1x,F5.3)') Iyr,qind
      ENDIF
C
C     Calculate the harvest guideline
      HG = 0
      Cpred = 0
      Eff = 0
      ObsEff = 0
C
      HG = HarvGuide(Isim,Iyr)
      ObsEff = HG/Bhat
C     Set the minimum catch of 2000t for bait
      HG = MAX(HG,2000.0)
C      HG = MAX(HG,0.0)
C
C     Make sure the catch does not exceed the biomass
      IF(HG.GE.B1p(Isim,Iyr-1)) THEN
        HG = B1p(Isim,Iyr-1)*0.95
      ENDIF
C
C     Save the F values
      Eff = HG/B1p(Isim,Iyr-1)
      IF(OFLopt.EQ.1) THEN
       Eff = (Emsymax-Emsymin)/(Nsim-1)*(Isim-1) + Emsymin
      ENDIF
C      WRITE(*,'(F9.1,1x,F8.2,1X,F9.1,1x,F9.1)') Bhat, EpsB(Iyr-1),HG, 
C     +                                          B1p(Isim,Iyr-1)
C      WRITE(*,'(I3,1x,A4,1x,F8.1)') Iyr,"HG", HG
C
C     SENSITIVITIES
C     3. Multiple fleets - Multiply HG by distribution
      IF(MultFleet.EQ.1.0) THEN
C       Don't forget to still set the minimum to the bait catch
        HG = MAX(HG,0.0)
C        HG = MAX(HG,2000.0)
C       Multiply by DISTRIBUTION (0.87)
        HG = HG * 0.87
C       And to make sure the catch does not exceed the biomass
        IF(HG.GE.B1p(Isim,Iyr-1)) THEN
          HG = B1p(Isim,Iyr-1)*0.95
        ENDIF
C      WRITE(*,'(I3,1x,A4,1x,F8.1)') Iyr,"HG2", HG
      ENDIF
C	  
C     Start calculation of mean age of the catch
      AgeC(Isim,Iyr) = 0
C     Solve for the catch
      IF (HG.LE.0) THEN
       FF = 0
C       WRITE(*,*) "Yes"
      ELSE
       Fmin = 0
       Fmax = 2.996
       FF = 0
C       FF = -1*LOG(1-Eff)
C      Use a bisection algorithm to guarantee that the HG is taken
       DO 11000 II = 1,23
        FF = (Fmin+Fmax)/2.0
        CPred = 0
        ACcount1 = 0
        ACcount2 = 0
        DO 11100 Age = 0,Nage
         Z = Selex(Age)*FF+NatM
         Cnum = N(Iyr-1,Age)*Selex(Age)*FF/Z*(1-exp(-Z))
         Cpred = Cpred  + Weight(Age)*Cnum
         ACcount1 = ACcount1 + Cnum*Age
         ACcount2 = ACcount2 + Cnum
11100   CONTINUE 
        IF (Cpred.GT.HG) THEN
         Fmax = FF
        ELSE
         Fmin = FF
        ENDIF
11000  CONTINUE        
      ENDIF 
C      WRITE(*,'(A4,1x,F5.3)') "FF", FF
      F(Isim,Iyr-1) = Eff
      AgeC(Isim,Iyr-1) = ACcount1/ACcount2
      IF(Cpred.LE.0.0001) THEN
        AgeC(Isim,Iyr-1) = 0.000
        F(Isim,Iyr-1) = 0.000
      ENDIF
C      WRITE(*,*) "F", F(Isim,Iyr), Isim, Iyr
      Catch(Isim,Iyr-1) = Cpred
C
C     Update the dynamics
      ANcount1 = 0      
      ANcount2 = 0
      AgeN(Isim,Iyr) = 0
      N(Iyr,Nage) = N(Iyr-1,Nage)*exp(-Selex(Nage)*FF-NatM)
      Z = Selex(Nage-1)*FF+NatM
      N(Iyr,Nage) = N(Iyr,Nage) + N(Iyr-1,Nage-1)*exp(-Z)
      DO 20000 Age = 1,Nage-1
       Z = Selex(Age-1)*FF+NatM
       N(Iyr,Age) = N(Iyr-1,Age-1)*exp(-Z)
       ANcount1 = ANcount1 + N(Iyr,Age)*Age
       ANcount2 = ANcount2 + N(Iyr,Age)
C       WRITE(*,'(I4,1x,I2,1x,F5.3,1x,F5.3)') Iyr, Age, FF, Z
20000 CONTINUE
C
C     SENSITIVITIES
C     3. Multiple fleet - Update dynamics
      IF(MultFleet.EQ.1.0) THEN
        ANcount1 = 0      
        ANcount2 = 0
        AgeN(Isim,Iyr) = 0
        N(Iyr,Nage) = N(Iyr-1,Nage)*exp(-Selex(Nage)*FF-NatM)
        Z = Selex(Nage-1)*FF + NatM + 
     +      SelxMC2(Nage-1)*0.2 + SelxPNW(Nage-1)*0.1
C        Z = MIN(Z,3.4)
        N(Iyr,Nage) = N(Iyr,Nage) + N(Iyr-1,Nage-1)*exp(-Z)
        DO 55000 Age = 1,Nage-1
         Z = Selex(Age-1)*FF + NatM + 
     +       SelxMC2(Age-1)*0.2 + SelxPNW(Age-1)*0.1
C         Z = MIN(Z,3.4)
         N(Iyr,Age) = N(Iyr-1,Age-1)*exp(-Z)
         ANcount1 = ANcount1 + N(Iyr,Age)*Age
         ANcount2 = ANcount2 + N(Iyr,Age)
C         IF (Iyr.GT.530) THEN
C           WRITE(*,'(I4,1x,I2,1x,F5.3,1x,F5.3)') Iyr, Age, FF, Z
C         ENDIF
55000   CONTINUE
      ENDIF      
C
C     Calculate mean age of the population
      AgeN(Isim,Iyr) = ANcount1/ANcount2
      IF(N(Iyr,Age).LE.0.0001) THEN
        AgeN(Isim,Iyr-1) = 0.000
      ENDIF	  
      IF(Iyr.EQ.Nyear) THEN
       IF(N(Iyr,Age).LE.0.0001) THEN
        AgeN(Isim,Iyr) = 0.000
       ENDIF	 
      ENDIF
C
C     Compute the SSB and 1+ biomass
      SSBV = 0
      Bio1P = 0
      DO 21000 Age = 1,Nage
C       SSBV = SSBV + Weight(Age)*N(Iyr,Age)
       Bio1P = Bio1P + Weight(Age)*N(Iyr,Age)
21000 CONTINUE
C
      DO 22000 Age = 2,Nage
       SSBV = SSBV + Weight(Age)*N(Iyr,Age)
22000 CONTINUE
C
      SSB(Isim,Iyr) = SSBV
      B1p(Isim,Iyr) = Bio1P
C
C     Get the recruitment ((EQUATION 2))
C
C     Base case is log(R/S) model ,
      RecR = SSBV/1000 * EXP(AlogRPS + BlogRPS*SSBV/1000 + 
     +       PhilogRPS*V(Iyr))
C     the log(R) model, dome-shaped      
      IF(Rtype.EQ.1) THEN
       RecR = EXP(Alpha + Beta1*SSBV/1000 + Beta2*(SSBV/1000)**2 +
     +        Phi*V(Iyr))
      ENDIF
C     The log(R) model, flat-top
C     First, set the upper limit of the recruitment curve
      SSBflat = -Beta1/(2*Beta2)
      IF(Rtype.EQ.2) THEN
       IF((SSBV/1000).GT.SSBflat) THEN
        RecR = EXP(Alpha+Beta1*SSBflat+Beta2*(SSBflat)**2+
     +         Phi*V(Iyr))
       ENDIF
      ENDIF
C     A simple Ricker model with no temperature term
      IF(Rtype.EQ.3) THEN
       RecR = SSBV*EXP(RickerA + RickerB * SSBV)
      ENDIF	   
      IF(Rtype.EQ.4) THEN
       RecR = SSBV*EXP(RickerA + RickerB * SSBV)
      ENDIF	 
      RecR = RecR * EXP(Eps(Iyr)-(SigmaR**2)/2.0)*1000
      N(Iyr,0) = RecR
      RecRep(Isim,Iyr) = RecR
C      WRITE(*,*) Isim,Iyr,SSBV,Recr
C      
      RETURN
600   FORMAT(I4,1x,100(F12.3,1x))      
      END
C
C ===============================================================================================
C
      SUBROUTINE Summary()
C
      IMPLICIT NONE
      INCLUDE "SARDINE.INC"
C      
C     Local variables
      INTEGER Isim,Year
C
      DO 10000 Year = 0,Nyear
       WRITE(14,600) Year,(B1p(Isim,Year),SSB(Isim,Year),
     +                     RecRep(Isim,Year),Catch(Isim,Year),
     +                     F(Isim,Year),ObsF(Isim,Year),Vrep(Isim,Year),
     +                     AgeN(Isim,Year),AgeC(Isim,Year),Isim= 1,Nsim)
10000 CONTINUE     
C
C     These are only for debugging. Comment them once the model works      
C      WRITE(14,*) F(1,2), F(1,3), F(1,4)
C      WRITE(*,*) F(1,2), F(1,3), F(1,4)
C      WRITE(14,*) Ninit
C      WRITE(14,*) N(1,0), N(1,1), N(1,2), N(1,3), N(1,4), N(1,5)
C      WRITE(14,*) SSB(1,0)
C      WRITE(14,*) V(1),V(2),V(3),V(4),V(5),V(6),V(7),V(8),V(9),V(10)
C
C      DO 11000 Year=0,Nyear
C       WRITE(14,*) Eps(Year)
C11000 CONTINUE
C	   
      RETURN
600   FORMAT(1x,I5,1x,100(F13.2,1x,F12.2,1x,F14.2,1x,F12.2,1x,F6.4,1x,
     +       F6.4,1x,F7.4,1x,F6.4,1x,F6.4))
    
      END
C
C ===============================================================================================
C
      REAL*8 FUNCTION HarvGuide(Isim,Iyr)
C
      IMPLICIT NONE
      INCLUDE "SARDINE.INC"
      
      REAL*8 HGt, Emsy,IEmsy,Frac,Imovave
      INTEGER Isim, Iyr
C
C     Calculate moving average of SST
      Imovave = (I(Iyr-1)+I(Iyr-2)+I(Iyr-3))/3
      IF(Iyr.eq.1) THEN
        Imovave = I(Iyr-1)
      ENDIF
      IF(Iyr.eq.2) THEN
        Imovave = (I(Iyr-1)+I(Iyr-2))/2
      ENDIF
C      WRITE(*,*) Imovave
C
C     Calculate Emsy
C     By default (OFLopt=0), set Emsy=Emsymin
      Emsy = Emsymin 
      IF(OFLopt.EQ.1) THEN
       Emsy = (Emsymax-Emsymin)/(Nsim-1)*(Isim-1) + Emsymin
      ENDIF
C
C     Calculate the OFL
      OFL = Emsy*Bhat
      IF (OFLopt.EQ.2) THEN
       Emsy = 15.28001 - 2.80651*Imovave + 0.16232*Imovave**2 - 
     +        0.0028564*Imovave**3
       IEmsy = Emsy
       Emsy = MAX(Emsy,0.0)
C       WRITE(*,'(A6,1x,F6.4)') "Emsy1",Emsy
       Emsy = MIN(Emsy,0.241)
C       WRITE(*,'(A6,F6.4)') "Emsy2",Emsy
       OFL = Emsy*Bhat	  
      ENDIF
C      WRITE(*,'(A6,1x,F6.4,1x,F10.5)') "Emsy2",Emsy,I(Iyr-1)	  
C
C     Calculate the HG
C     By default (HCRopt=1), set the HG = Emsy
      HGt = OFL
      ObsF(Isim,Iyr-1) = Emsy
C     It can also be set with a cutoff and maxcatch and pre-specified Fraction=FracMin
C
C     Cutoff can be looped over 
      IF (CutOpt.EQ.1) THEN
       Cutoff = 50000*(Isim-1)
      ENDIF
      IF (HCRopt.EQ.2) THEN
       Frac = FracMin
       HGt = Frac*(Bhat-Cutoff)
       HGt = MIN(HGt,Maxcatch)
       HGt = MIN(HGt,OFL)
       HGt = MAX(HGt,0.0)
       ObsF(Isim,Iyr-1) = Frac
C       WRITE(*,'(F10.1,1x,F10.1,1x,F10.1)') Bhat,HGt,OFL
C     Or include a temperature term to determine Fraction
       ELSE 
        IF (HCRopt.EQ.3) THEN
         Frac = Emsy
         Frac = MIN(Frac,FracMax)
         Frac = MAX(Frac,FracMin)
         HGt = Frac*(Bhat-Cutoff)
         HGt = MIN(HGt,Maxcatch)
         HGt = MIN(HGt,OFL)
         HGt = MAX(HGt,0.0)
         ObsF(Isim,Iyr-1) = Frac
        ENDIF
      ENDIF
C	  
      HarvGuide = HGt
C
C      WRITE(*,'(A4,1x,F6.4,1x,A4,1x,F9.1)') "HGt", HGt, "OFL", OFL
      RETURN
      END
C
C ===============================================================================================
C
      INCLUDE "COMMON.FOR"